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SUMMARY 

The use of the interaction integral to compute stress intensity factors around a crack tip requires 
selecting an auxiliary field and a material variation field. We formulate a family of these fields 
accounting for the curvilinear nature of cracks that, in conjunction with a discrete formulation of 
the interaction integral, yield optimally convergent stress intensity factors. We formulate three pairs 
of auxiliary and material variation fields chosen to yield a simple expression of the interaction integral 
for different classes of problems. The formulation accounts for crack face tractions and body forces. 
Distinct features of the fields are their ease of construction and implementation. The resulting stress 
intensity factors are observed converging at a rate that doubles the one of the stress field. We provide 
a sketch of the theoretical justification for the observed convergence rates, and discuss issues such 
as quadratures and domain approximations needed to attain such convergent behavior. Through 
two representative examples, a circular arc crack and a loaded power function crack, we illustrate 
the convergence rates of the computed stress intensity factors. The numerical results also show the 
independence of the method on the size of the domain of integration. 
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1. INTRODUCTION 

The stress field near the tip of a loaded crack is singular under the assumption of linear elastic 
fracture mechanics. The coefficients of the asymptotic stress field, known as the stress intensity 
factors, play a key role in characterizing the magnitude of the load applied to the crack and 
predicting its propagation. 

Given the stress singularity and the poor accuracy in pointwise evaluation of the stress field, 
it is often impossible to extract the stress intensity factors directly from numerical solutions. 
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unless a higher-order method to compute the elastic field is adopted, such as those proposed 
in Liu et al. [1], Shen and Lew [2], and Chiaramonte et al. [3]. 

As a result, path and domain integral methods to extract the stress intensity factors have 
been created precisely to circumvent this limitation. A method of this kind typically formulates 
the expression of the stress intensity factors as functionals of the solution, thus enjoying a 
higher order of convergence than the one of the elastic field itself. Predominant methods of 
this kind have been constructed based on the J-integral [4] and the interaction integral [5]. 

In the context of linear elastic fracture mechanics, the J-integral is identified with the 
system’s elastic energy release rate; the elastic energy that would be released per unit length of 
crack extension in the tangential direction. This integral and related ones for the computation 
of fracture-mechanics-related quantities have been elaborated by Eshelby [6], Rice [4], Freund 
[7], and many others. Shih et al. [8] derived the expression for the energy release rate of 
a thermally stressed body in the presence of crack face traction and body force. A general 
treatment of such conservation integrals, including those expressed in a form of domain 
integrals, can be found in Moran and Shih [9, 10]. The domain form is better suited and 
more accurate for numerical computation. Nevertheless, in the case of mixed-mode loading, J 
is a quadratic function of all three stress intensity factors [11]; therefore, additional integrals 
are needed to determine the three quantities individually, e.g., as done by Chang and Wu [12] 
for non-planar curved cracks. 

In contrast, the interaction integral, or the interaction energy integral, is able to yield the 
three stress intensity factors separately. This method is based on the J-integral by superposing 
the elastic field of the loaded body and an auxiliary field with known stress intensity factors. 
The auxiliary field does not need to satisfy the elasticity equations but must resemble the 
asymptotic solution of a cracked elastic body corresponding to one of the three loading modes 
(e.g., plane-strain mode I or mode II, or anti-plane mode III). Therefore the auxiliary field, for 
straight cracks, is normally chosen to be the asymptotic solution, as found in [13, 14]. Doing 
so readily yields the stress intensity factor of the actual field for the chosen mode. Along with 
the auxiliary field, the interaction integral requires the construction of a vector field, named 
the material variation field, which indicates the velocity (variation) of points in the reference 
configuration as it is deformed into a domain with a longer crack. Under mild conditions, the 
value of the interaction integral does not depend on this choice, but a good choice of material 
variation can simplify computations. For example, the interaction integral is computed by 
integrating over the support of the material variation field, so it is convenient to choose material 
variation fields with small and compact support. While developing auxiliary and material 
variation fields for straight cracks (planar cracks in three dimension) is an amenable task, 
doing so for curvilinear cracks (non-planar cracks in three dimension) poses several challenges. 
In the following paragraphs we provide a short review of the effort related to the computation 
of stress intensity factors with the use of the interaction integral. 

Earlier methods to compute the stress intensity factors with the interaction integral involved 
path integrals, such as those in Stern et al. [5] and Yau et al. [15]. The method was later 
generalized to a straight-front crack in three-dimensions by Nakamura and Parks [16] and 
Nakamura [17]. A curved crack front introduces additional terms, since the popular plane- 
strain modes I and II auxiliary fields no longer satisfy the compatibility and the equilibrium 
conditions. These additional terms were accounted for in Nahta and Moran [18] for the 
axisymmetric case, and in Gosz et al. [19] for general planar cracks in three dimensions. Kim et 
al. [20] adopted auxiliary fields corresponding to penny-shaped cracks, as well as conventional 
plane-strain and anti-plane ones. An alternative approach was given by Daimon and Okada 
[21] who adopted a compatible auxiliary field and accounted for its lack of equilibrium by 
superposing a numerically computed displacement field with the finite element method. A 
study of the effect of omitting some terms accounting for the curved front is given by Walters 
et al. [22]. 

More recently the method of [19] was adapted for non-planar cracks in Gosz and Moran [23], 
the latter of which is arguably a milestone in the development of domain integral methods to 
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extract stress intensity factors from curvilinear cracks in 2D and non-planar cracks in 3D. 
In [23], the method constructs the auxiliary fields through the use of curvilinear coordinates 
and their corresponding covariant basis to account for the crack curvature. This procedure 
constructs the auxiliary fields by juxtaposing the components of the stress fields of [13] with the 
described basis. In [24], Sukumar et al. implemented the method of [23] in combination with an 
extended finite element method in a three-dimensional setting [25] and a fast-marching method. 
The main drawback of the curvilinear coordinates of [23] is the need to perform boundary 
integrals over both the real crack surfaces and a pair of fictitious crack surfaces. In [26], 
Gonzalez-Albuixech et al. proposed another curvilinear coordinate system that can eliminate 
the integration on the fictitious crack surfaces, and in this way facilitate the computation. 

In [27, 26], Gonzalez-Albuixech et al. studied the properties of the aforementioned methods 
in two- and three-dimensions, respectively, but not with all terms arising from the derivation 
of the interaction integral. Such omission of terms may have contributed to the observed slow 
convergence and occasional divergence. This observation confirms the statement of [23] that 
all terms arising from the lack of compatibility and equilibrium of the auxiliary field have to 
be taken into account. 

The domain version of the interaction integral has also been generalized to functionally 
graded materials [28, 29]. 

In this paper we present a suite of auxiliary fields and material variations fields. By pairing 
two constructs of material variation fields and two constructs of auxiliary fields, we create 
two kinds of interaction integral suitable for curvilinear cracks and for situations in which 
body forces and crack face tractions are present. One kind of interaction integral is suited 
for applications where crack faces are loaded (e.g. hydraulic fracturing), and the other one is 
best suited for applications where body forces are non-zero (e.g. thermally loaded materials). 
Moreover, no fictitious crack face is needed, a major simplification to the predominant method 
in the literature. 

One of the two choices of the material variation fields has a constant direction pointing to 
the direction of the crack growth, a straightforward choice adopted by most authors, and the 
other has a direction that is tangential to the crack near the crack tip, similar to that proposed 
in [25]. For a curved crack, this second choice necessarily coincides with the first one only at 
the crack tip. A key advantage of the material variation fields we introduce here is their ease 
of construction which is reflected in their straightforward implementation in computer codes. 
Moreover, in contrast to many of the existing constructions, the magnitude of the material 
variation fields is mesh independent. This mesh independence contributes to the observed 
optimal rate of convergence. 

The two auxiliary fields are constructed from the well-known asymptotic solutions of a 
straight crack. Both fields respect the discontinuity introduced by the crack, thus avoiding 
the evaluation of integrals over fictitious crack faces as the method in [23] does. One of the 
auxiliary fields is obtained by “extending” the asymptotic solutions past the range [— 7r,7r], 
and hence satisfies equilibrium, compatibility, and the constitutive relation. The resulting 
interaction integral expression then yields a term on the crack faces, even in the absence of 
crack face traction. The second auxiliary field is an incompatible strain field. It is obtained 
by first mapping a straight crack to the curved crack near the crack tip, and using this map 
to push forward the strain field of the straight-crack asymptotic solution. Then, by suitably 
rotating the strain tensor at each point, we obtain an auxiliary strain field that is traction-free 
at the curved crack faces. This is useful for problems in which crack faces are traction-free. 
In fact, if this auxiliary field is used in combination with the tangential material variation 
field, the crack face integral vanishes, resulting in a significantly simplified expression for the 
interaction integral. 

We showcase the convergence of the stress intensity factors obtained with the proposed fields 
for a set of representative examples computed with two different finite element methods. In 
all cases, the stress intensity factors converge with a rate that doubles the rate of convergence 
of the strains. We also numerically demonstrate the independence of the computed stress 
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intensity factors from the chosen support for the material variation field. Although the 
numerical examples adopt finite element methods to obtain an approximate solution to the 
elasticity problem, the numerical implementation of the interaction integral with the new 
fields is general, and can be used in conjunction with any numerical method for the solution 
of the governing equations (e.g. finite difference, finite volume, boundary integral equations, 
isogeometric analysis, and meshless methods). 

The paper is organized as follows. We first state the problem that we seek to solve in §2. We 
then proceed in §3 to present the interaction integral with the description of the new material 
variation and auxiliary fields. In the same section we justify that the proposed forms of the 
interaction integral are well-defined. A numerical approximation of the interaction integral is 
presented in §4 with remarks on its expected convergence. The last part of §4 provides a step- 
by-step recapitulation of the method suited for the reader interested in a concise presentation. 
In §5 we verify the computation of the stress intensity factors against analytical solutions 
for two problems: a circular arc crack and a power function crack. Throughout the paper 
we included sections titled ^^Justification'’’ which contain sketches of proofs for some of the 
assertions we make, and they are not essential for the description of the methods in this 
paper. 


2. PROBLEM STATEMENT 

We present next the problem statement which consists of the evaluation of the stress intensity 
factors following the solution of the elasticity fields for a cracked solid. 

2.1. Elasticity Problem 

We consider a body C undergoing a deformation defined by the displacement field u. 
We assume B to be an open and connected domain with a (piecewise) smooth boundary dB. 
We represent the crack with a twice differentiable, simple and rectifiable curve ^ C B and 
denote its crack faces with . The cracked domain is given by = B\'t^. The boundary 
of B<^ is the union of the crack faces and the boundary of B, namely dB<^ = dBL}'l^±. Let 
dB’^ be decomposed into drB’^ and dd.B’^ such that drB<^ 3 drB<^ L) ddB<,f = dB<^, and 
drB<^ n ddB<^ = 0. Tractions t and displacements u are prescribed over drB<^ and ddB<^, 
respectively, while a body force field b is applied over B<^. Let Xt denote any one of the 
two crack tips. We denote by n the unit external normal to B, as well as the unit external 
normal to each one of the two faces of the crack. Figure 1 shows the schematic of the problem 
configuration. 


t 
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We confine ourselves to planar linear elasticity in the context of infinitesimal deformations. 
The elasticity problem statement reads: Given b , u, and t find u : B^if —> such that 


V • cr(VM) + 6 = 0, (la) 

u = u, on ddB<^g, (lb) 

a-{\7u)n = t, on drB<^, (Ic) 

where cr is the stress tensor. This is given by 

cr(VM) = C : Vm, (2) 


and 


for plane strain, 
for plane stress. 


. A, 

C = A 1 G 1 + \ = ^ 2X^1 

i A + 2/i 

The constants A and /i are Lame’s first and second parameters, respectively, 1 is the identity 
second-order tensor, and I is the fourth-order symmetric identity operator given by 


H — ^ “t“ G G (^k G e/, 


where { 61 , 62 } is a Cartesian basis, and an index repeated twice indicates sum from 1 to 2 in 
such index. 


2.2. Crack Tip Coordinates and Stress Intensity Factors 

To aid the definition of the stress intensity factors we first introduce a system of coordinates 
and a family of vector bases. 



Figure 2. Description of local basis and coordinates 

Let (r,!)) be crack tip polar coordinates as shown in Fig. 2, and let Bp{x) = 
{a; e I |a; — a:| < p| be the open ball of radius p > 0 centered at x. The radial coordinate 
is defined as r{x) := \x — Xt\. Let T ■. r ^ {x & Io\\x — Xt\ = r} he & description of part of 
the crack parametrized by the distance to the crack tip. We set the domain of T to be 
[0,p] with p > 0 such that Bp{xt) C B and that F' 7 ^ 0 over its domain of definition. As a 
consequence, F is bijective for r G [0,p]. We re-iterate that the crack is assumed to be a 
twice continuously differentiable curve such that F G C'^(M(}';M^), where M)} = {0} UM+. By 
convention, the possible values of the (r, i9) coordinates for points in Bp{xt) that we will use 
belong to 

Dp = {ir,'d) X M I -TT - C(r) < -d < TT - C(i’)}, 

where C(i’) is the angle between the vector F(r) — Xt and r'(0). In other words, C(r) is the 
angle subdued by (a) the tangent at the crack tip and (b) the secant line passing through the 
crack tip and r(r). Figure 3 shows the values of the coordinate d for the particular case of a 
circular arc crack. Lastly let gi(r),r G [0,p] be the right-handed orthonormal bases induced 
from the mapping F such that gi{r) = —F'(r)/|r'(r)|, see Fig. 2. 
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—TT + max ( i9 TT + max ( 


Figure 3. The i9-coordinate with the branch cut along the crack. 


Throughout this manuscript, we assume that there exist unique real numbers Kj, Kjj such 
that 

u = Kju^ + Kiiu^^ + 1x5, Us € , (3) 

where Kj and Kjj are the modes I and II stress intensity factors, respectively, and G 

II^(B<jf;M^) \ II^(5^;M^) are the asymptotic displacement solutions of these modes. 

The problem of evaluating the stress intensity factors is: Given a solution u of Problem (1), 
compute K/ and Kjj as dehned in (3). 

Remark (Explicit evaluation of the behavior at the crack tip). Whenever /3 = Vtt € 
the stress intensity faetors can be alternatively defined as: 

Kj = lim : ^2 G) §2, 

(4) 

Kjj = hm V27rrcr(/3)|^^o : Si <8 92 - 

r —>0 

The calculation of the stress intensity factors by evaluating the limits in (4) with a numerical 
solution often leads to poor results. In fact, few methods are eapable of aeeurately resolving 
the singularity in the stress field. Therefore the predominant methods to eompute the stress 
intensity faetors are based on the approximation of the interaction integral, a formulation that 
avoids pointwise evaluation of the stress field in the region with the singularity. We proceed to 
introduce the interaetion integral in §5. 


3. INTERACTION INTEGRAL 

In the sequel we first define the interaction integral functional between any two admissible 
fields alongside a concise justification of this definition (see §3.1 ). In §3.2 we specialize it to 
the case in which one of the fields is the solution to the elasticity problem of interest. This 
specialization results in a formulation that possesses the following properties: (a) it does not 
involve second derivatives of the discrete approximation to the exact displacement field; (b) it 
is a problem-dependent functional because it uses the prescribed tractions and body forces; and 
(c) it can be further simplihed, depending on the problem, by carefully choosing the so-called 
material variation and auxiliary helds, to be described later. We perform the simplihcation 
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in §3.3 and §3.4, and obtain three different formulas. In §3.5 we provide guidelines on which 
specific formula for the interaction integral to choose for a given application. 

3.1. Definition of the Interaction Integral Functional 

The interaction integral involves two elasticity fields, one with known stress intensity factors, 
and the other whose stress intensity factors we are interested in evaluating. 

Let us introduce the interaction energy momentum tensor S : x defined 

as 

S {f3\(3^) = w (/3“, ^3") 1-(3^^a- {(3^) - (3^ 
where w : x —)• M is given as 

W (/3“, (3^) = ^ [<t(/3“) :/3^ + a (/j") : /3“] . (5) 

Assuming the same constitutive relation (2) for both fields, (5) simplifies to 

w(/3^,/3^) =ct(j 3-):/3^ = a (j3^) :/3F 
Additionally, let the set of material variations be dehned as 

M = {(I7 e |(57 = Bp{xt), 6 j{xt) = 91(0)} . (6) 

Finally, let 

= span{VM^ {B^;R‘^^‘^) , 

and for any tensor field (3 G define Ki[P'] and Kjj[(3] such that 

(3 = Ki[f3]Vu^ + Kii[f3]Vu^^ + f3s (7) 

with /3g G In particular, if /3 = Vtt, for u being a solution of Problem (1), 

then K[[f3] and Kii[f3] are the stress intensity factors of u. However, Kj[l3] and Kii[f3] are 
still defined for any /3 G that is not the gradient of a displacement field. For convenience, 
regardless of whether (3 is or is not the gradient of a displacement field, we will refer to Ki[l3] 
and Kii[0\ as the stress intensity factors of /3. 

We define the interaction integral functional X : x x AI —?• M as 


X [/3“, (3\ 6j] = [ (57 • S (/3“, f3^) n dS 


' Bp{xt)Y^ 


S (/3“, /J") : V(57 + V • S (/3“, (3^) • < 57 ] dV, 


where, for convenience, we let denote '^±1^ Bp{xt). The value X [/3“,/3^, (57] is the 
interaction integral between /3“ and (3^. The relation between the interaction integral and 
the stress intensity factors of [3°', [3^ G is 

X[(3\l3\5-i] =g{Ki[^'^]Ki [/3'] + [/j"]) (9) 

for any (57 G A4, where rj is a material constant defined as 


^ A © 2/r 
2/i(A + /i) 

rj = < 

2(A + /i) 

. fi(3X © 2/i) 


for plane strain. 


for plane stress. 


It follows from (9) that, if we are interested in finding Kj[l3] (or Kii[(3]), we must generate 
an auxiliary tensor field (or f3ff^) G satisfying Ki[f3Y^^] = 1, Kii[f3Y'^] = 0 (or 


Prepared using nmeauth.cis 



M. M. CHIARAMONTE, Y. SEEN, L. M. KEER, AND A. J. LEW 


= 0) = 1). In this case, (9) implies that 


Kiji[f3] 


V 


Notice that the interaction integral can be regarded as a tool to extract the singular parts 
of fields as it follows from (8). The regular part /3g of either field (cf. (7)) does not 

contribute to the value of the interaction integral. 

Justification (Equation (9)). Consider 13°', [3^ G and r > 0. Notice that the two terms 
in the volume integral of (8) form an exact divergence. Applying the divergence theorem 
on {Bp{xt) \ ^) \ Br{xt) reveals that the integration in (8) over {Bp{xt) \ ^) \ Br{xt) and 
\ Br{xt) add up to an integral over dBfixt). It then follows that 


I [13°, f3\ , 57 ] = lim / 57 • S {f3°,f3^) n dS, (10) 

with n here is also used to denote the outward unit normal to dBr{xt), since the rest of the 
terms vanish as r —> 0. 

To proceed in showing that (10) implies (9), we write (3° = f3^ + (3°g and f3^ = (3^ + l3^g, 
where f3^,f3^ € spanjVM^, and f3g,f3g € It is straightforward to show 

that 

J [f3°r,(3^T, Si] = lim [ <57 • S {(3°r, /3^) n dS = rj {Kj[(3°]Kj [f3^] + Kii[^°]Kii ) . 

Therefore, it remains to show that i[f3g, /3^, ^ 7 ] = 'i[l3'^, f3g, < 57 ] = i[l3°g, f3g, < 57 ] = 0. To this 
end, we first define 


ir[(3°,(3\Si]:= [ Si-i:{f3°,f3’>) 

J dB-r-ixA 


n dS. 


Then we invoke the Cauchy-Schwarz inequality and the explicit expression of /3^ to obtain 


ir [/3s) ^T> <^ 7 ] I < C\\l3°g\\g2yQB^(xt)] < C'll/3sllL2[9B,.(xt)]) 

[f3°g,(3%,5i] I < C||/3^|U2[aB.(x71 ||/3s|L2[aB.(x7] > 


Ir 


( 11 ) 


where C > 0 is independent of r. 

To continue, we need to invoke a trace inequality with a scaling of r for any / G iJ^[i?p(a:i) \ 
and r G (0,p), 

ll/llL79B,(xt)] < C'?'^^^ll/llRi[Bp(xt)\'5f], (12) 

where C > 0 is independent of / and rfi 


t To prove (12), we first write 


/:= 


Jb. 


fdn, / := / - /. 


Then with a form of Poincare’s inequality and a scaling argument, 


IrlL^as.CxD] “ 

On the other hand, since \\f\\L2[gB^(xt)] = (7^27rr)72 , we have 

l|7lL2[aB7x71 = ll7lL7np(^71 ^ Il7llr7n.(-dl ■ 

Adding these two inequalities yields (12). 
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With (12), we then proceed to simplify (11): 


Thus, as r —)■ 0, both %■ [^ 5 ,/ 3 ^,( 57 ] and Xr [/^s,/S tj'^t] tend to 
the first two slots of X^, we also have limr_).oX[/3^,/ 35 , ^ 7 ] = 0 . 


zero. From symmetry in 

□ 


Remark (Relation to the energy release rate). The interaetion integral functional is directly 
related to the energy release rate Q : x A1 —> M [SO], which can he defined as 

g[[3,6j] = lim / 6j ■ i:{f3)n dS, (13) 

where S : —>• is Eshelby’s energy momentum tensor [31] and n is used to denote 

the outward unit normal to dBr{xt). For linear elastic materials, Eshelby’s energy momentum 
tensor takes the form 

S(/3) = ^^(/3):/3 1-/3T^(/3). 

The above is related to the interaction energy momentum tensor by the following relation 

S(/3“,^^) =S(/3“ + ^^)-S(/3“)-S(/3*') ^ {f3^) = {f3\ f3^) . 

Comparing (13) and (10) and exploiting the linearity of the constitutive relation, we have the 
following relation between the interaction integral functional and the energy release rate 

X [/3“, /3^ 57] = g [r + ( 3 \ 6 j] - g[[ 3 \ 57] - e [/3', ^7] • (14) 

After replacing with (7) and evaluating, the limit in (13) gives the widely known relation 

g[^,5-t]=p{Ki[(3f + KuW). 

We can alternatively recover (9) by replacing this relation into (14). Thus, (9) can be justified 
by the direct evaluation of the limit in (10) or by its relation to the energy release rate. R is 
worth noting that while g is a non-linear functional in (3, X is linear in both /3“ and /3^. 

Remark (Constraints on A4). In (13), Sj is understood as a variation of material points that 
represents unit crack advancement. Thus, in order to compute the energy release rate, we must 
enforce Sj = QiiO) at Xf. This constraint is the justification behind the definition of M. in (6). 

3.2. Problem-dependent Interaction Integral Functional 

We present here a functional X that takes the same value as X when /3“ coincides with 
the solution of Problem (1). Namely, if /3“ = Vtt where u satisfies (la) and (Ic), then 
X[f3°', (3^, (^ 7 ] = X[(3°‘,I3'’, 57 ] for any f3^. Therefore, in this case X[/3“, (3^, ( 57 ] is the interaction 
integral between /3“ and (3^. 

The motivation behind introducing X is to formulate a functional defined over gradients 
of displacement fields that belong to classical finite elements spaces, namely, a functional for 
which no second derivatives of a numerical solution are needed. This is possible because, in 
contrast to X, X does not involve derivatives of /3“. 

Expanding the divergence in (8) and substituting with (la) and (Ic) for /3“ = Vit yields 

l[[3\(3\5-i]= [ 57-t(/3“,/3') 

- [ ■.VSj + Sj-X{[3^,(3^)] dV, 
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where 


and 


T (/3“,/3'’) = w (/3“, (3^) n-l3^^a (/3*') n - /3'’ 

A (/3“, /3^) =13^ :Va- (/3'') - cr (^3“) : (V/3*' )^ - /3“ • cr (/3'’) + /J** . 


(15) 

(16) 


Remark (Indicial expression of relevant quantities). For ease of implementation, we provide 
here the indicial representation of S, r, and A (making use of Einstein’s repeated indeces 
convention), namely 


S,, (/3“,/3'>) = ahfdl S., - 4“ 

T, {[3°', 13^) = wm - - Pjiij, 

(/ 3 “)/ 3 ^) = Pmn'^mn,i ~ '^kjPkiJ ~ Pki^kj,j + Pki^k, 


where are understood as a{(3°'’^)ij. 

Notice that for each pair (^f3°',f3^), T(^f3°',f3^) and X(^j3‘^,f3^) are functions over and 
Bp{xt) \ ^, respectively. 

To reflect the lower regularity needed for /3“, we first define a finite partition of as a set 
{Ti,..., Tpf} for some N gN such that Ti is open for any i, R n Tj = 0 for any i j, and 
(Jj Ti = An example of such partition is a finite element mesh for B<g. Then, we can set 


= ^(3 € if /3|Ti G F[^ for any Ti in some finite partition of 

0 span {Vm^, 


and define I: SS°‘ x x Ad —)■ M as 


3l 


l[(3‘^,(3\S-f]= 5-f-T {(3^,(3^) dS 


[ [•£{(3^(3'^) :WSj + 6j-X{(3\(3^)]dV. 

JBp(xt)\^ '' -^ ^ 


(17) 


l2 


I 3 


In this way, /3“ can have discontinuities across a finite number of interfaces, as in a finite 
element solutions, and still have well-defined values at the crack faces (which a function in Lf 
may not have). Note that I is linear in (3^ but affine in /3“, and therefore, not symmetric with 
respect to them. The way I will be used is by setting /3“ to be either the exact solution or a 
numerical approximation to it, and j3^ to be an auxiliary field. 


3.3. The Fields 

We next proceed to construct the material variation and auxiliary fields that will enable the 
extraction of the stress intensity factors for curvilinear crack geometries. 


3.3.1. Material variation fields. The objective of this section is to construct vector fields Sj 
that belong to the set of material variations Ad, see ( 6 ). We provide two constructs, but any 
(I 7 € Ad could be used. 

We start from a general form S^{xt + re^) = q{r)t{r) where q : IR(|" —)■ Mjj" and t : —)■ 

The function q{r) represents the magnitude of the material variation field and is constructed to 
have support within Bp{xt). The function t(r) embodies the direction of the material variation 
field and is taken to satisfy |t(r)| = 1, Vr G Mq . 

The magnitude and direction of the two material variation fields that we propose only depend 
on r. We will thus abuse notation writing 5^{r) in place of d^{xt + re^). 
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The scalar function q{r) G is defined as 


q{r) 


/ 

1 , 

< fir), 

. 0 ^ 


Hr < pi, 
ii Pi < r < p, 
otherwise, 


(18) 


with /(r) being a fifth order polynomial and pi = p/4. Note that to construct higher-order 
methods the regularity of q{r), and thus the polynomial order of f{r), will have to be suitably 
adapted. 

In the sequel we list the two material variation fields: 


(1) Unidirectional material variation fields. The first field is designed to be constant within a 
distance pi from the crack tip. The field is then constructed as 

<57™'(r)=g(r)9i(0). (19) 

This field satisfies 

= 0, for r < pj. (20) 

Figure 4a shows its stream traces alongside a circular arc crack. 

(2) Tangential material variation fields. The second field is designed to be tangential to the 
crack and is given by 

Sj'rAN^r)=qir)g,{r). ( 21 ) 

This field satisfies 

,57^^^ • n = 0 on (22) 

The stream traces of < 17 "'"^^, for the particular case of a circular arc crack, are shown in 
Figure 4b. 



(a) Constant within \x — xt\ ^ Pi 



Figure 4. Streamtraces for the material variation fields: (a) (b) 57 ^^^. 


Remark (Regularity of Sj). Note that since T G (^^(Mg ) and q G C^(M(j"), both and 

S^tan ^/jg continuity requirement of A4, namely, both are in 
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3.3.2. Auxiliary fields. As discussed in §3.1, the objective is to construct tensor fields 

G (23a) 

such that 

Kj[f3r] = l, Kjj[f3r]=0, 

Kj[f3^r] = o, = 

For a crack that is straight near the tip, namely, ‘rf n Bp{xt) is straight, a natural choice is 
the strain fields of the solutions to pure modes I and II loading [13]. In fact, these solutions, 
appropriately scaled, satisfy (23b) and the regularity requirement (23a). Furthermore, the 
stress field cr(/3j™j) is divergence free and the fields /J/”// are compatible, i.e., 

V • =0 in Bp{xt) \ (23c) 

3$: Bp{xt) \ such that = V^> in Bp{xt) \ (23d) 

as they are indeed derived from gradients of vector fields. Additionally the stress field is 
traction-free on the crack faces: 

cr{(3Y^fj)n = 0 on (23e) 

These features allow for significant simplifications of the interaction integral functional in (17). 

For curvilinear cracks, however, analytically obtaining auxiliary fields with the same features 
is not generally possible, since a field that satisfies all conditions (23) is the solution of Problem 
(1) in the neighborhood of Xt for the given curvilinear crack geometry Instead, we will 
construct auxiliary fields that, although sufficiently regular and satisfying (23b), may violate 
(23c), (23d), or (23e). Needless to say that doing so hinders the simplification of the interaction 
integral functional, as discussed in §3.4. 

In the following we discuss two constructs of the auxiliary fields that satisfy (23a) and 
(23b): (1) we present a compatible with divergence-free stress field cr(/3j™/), but for 
which (t(/3j“^j) is not traction-free on the crack faces, and (2) then we introduce a variant of 
that is incompatible and whose stress field is not divergence free, but its stress field is 
traction-free on the crack faces. 

(1) Divergence-free and compatible (DFC) fields. We first construct an auxiliary field which 
satisfies conditions (23a), (23b), (23c), and (23d). 

To this extent consider the displacement fields obtained for a straight crack in pure 
mode I, II loading given by g^{A) ® 9^(0), where, for completeness, the 

components are recapitulated in §A. The auxiliary fields are then taken as 

d), (24) 

where for each r the domain of definition of A is [—tt — C(r),7r — ('(?')] as introduced in 
§2.2, rather than [—7r,7r]. 

(2) Traction-free (TF) fields. We now construct auxiliary fields such that (23a), (23b), and 
(23e) are satisfied. 

Consider the mapping ip: Dp —>■ [—7r,7r] of the angular component of the polar coordinate 
system introduced in §2.2. This mapping is designed to take a value of ±7r on the crack 
faces and can be constructed as 


(p{r,A) = d 3- ({r). 

Values of p are plotted for a circular arc crack geometry in Figure 5. 
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We then construct /S/as 


Vm 


I,II 


: 9 ,( 0 ) 0 9 ,( 0 ) 




9tir)^gAr). (25) 


This auxiliary field is well defined for r G [0, p), where C(r) is also well defined. Its values for 
r > 9 do not participate in the interaction integral, because of the support of Sj, and hence 
are immaterial. Below we show that jSj/j € span |0 and 

hence that it can be extended to a function j3jj/ G span 0 M^^^) = 

The inspiration behind this construct is to transport from the straight crack faces, 

on which cr{Vu^'^^) is traction-free, to the faces of the curvilinear crack, rotating 
and hence criVu^’^^), precisely by the angle between g^r) and 9^(0). This is generally a 
incompatible field with non-divergence-free stresses but traction-free crack faces. 


Justification (Traction-free property). We begin by computing the stresses from the 
constitutive relation (2) on both sides of (25). Let then • ]g-i- x (—7r,7r) —)■ 

denote cr{Vu^’^^), which are precisely the stress fields of a straight crack ( see §A) parallel 
to the local crack tip basis vector 9i(0). These stress fields are traction free along these 
straight faces, so ’^\r,±.TT)g 2 A)) = 0. Then, on we have 




rb//, 


r,(9(r,d)) : 9*(0) 0 9,(0) 


9,(0 0 9 , ( 0 « 






9,(0)-cr’ (r,±7r)9 -(0) 


9,(002 = 0 , 


where we used that p = ±7r on . 


(26) 


□ 


Justification (Regularity of /3"^^). It is not a priori apparent that 0^^ G , but it does. 
To prove (3^^ G first note that G spanjVM^, It is then enough to show 

that j3g := /Sjy/ —/3/J/ € H^{Bp{xt)]R^^Aj hence that it can be extended to a 
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function fUg G To this end, we write 

/ 3 s = + ((?■)) : 9^(0) (g) 9^(0) - Vu^’^\r,'d) : 9,(r) (g) 9^ (r)] 9,(r) (g) 

(27) 

where C(?') for r > 0 satisfies that 

^ r'(o) r(o)-r(r) 

cosC(r) = s(r) := - — 


|F(0)| |r(0)-r(r)r 


and C(0) = 0 = limr-s.o+ C(^)- Hence, 


/■// ^ I <3 . . s'(r) 

C (r) = ±— arccoss(r) = =F— , 

dr Vl - siry 


= Tt 


1 


r(o) • r(r) r(o) • [r(r) - r(o)] r(r) • [r(r) - r(o)] 


1 - - 


|r(r)-r(o)| 

r'(o) r(o) - r(r) 


|r(r)-r(o)|3 


|F(o)| |r(o)-r(r)| 


- 1/2 


which is well defined and continuous for any r > 0. A tedious calculation shows that 
C'(0) = limr_,.o+ Cir) is also well-defined and given by 


|C'(0)| = lv/-s"(0)| = 2 |p( 0)|2 {|r'(0)nr"(0)p - [r'(o) • r"(o)]"} 


1/2 


< OO. 


Hence, there exists C > 0 such that for all r G [0,/?], 

|C'(r)|<C, (28a) 

and thus 

|C(r)|<Cr. (28b) 

Here and henceforth C indicates a positive constant independent of r G [0, p], whose value 
may change from line to line. 

Next, as shown in §A, VM'^’'^^(r, d) = r“^/^/'^’^^(0), where G From 

(27), we can write 


fds = r [/(^ + CW) : 9i(0) (g)9j(0) - /(d) : 9 ,(r) (g) g-^/r)] 9 ,(r) (g) g^/r). 


where we have omitted the superscript I, II of /, as we shall do hereafter. It is then 
straightforward to show that /3g G 

It remains to show that 


df3s I dps 

dr ' r dD 


GL2 . 


We first examine 

r ^ [/'(^^ + C(0) : 9*(0)(g)gj (0) - /'(d) : 9 ,(r) (g) g^/r)] g,(r) (g) g^-(r) 

'iJ 


= E{ + CW) - ^ ® 9,(0) 

'ij 

+ /'(d) : [g, (0) (g) g^- (0) - g, (r) (g) (r)] jg^(r) (g) g^- (r). 

(29) 
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Since / is C°°, we apply (28b) and write 


ll/'(i^ + C(r-))-/'(i^)IU < \\f\\w^.^\ar)\<Cr. 


(30) 


On the other hand, we note that g;^(r) = —r'(r)/|r'(r)| and g 2 {r) = u) ■ g^{r) where 
w := —ei (g) 62 + 62 0 ei, differentiating with respect to r yields 


g'lif) 


r"(r) 

|r'(r)| 


r'(r) • r"(r) 
|r'(r)|3 


r'(r), 


g'^ir) = uj ■ g[{r), 


for r e [0,p]. Thus g'i{r) and 92(^) bounded, and hence 


Il9*(0) 0 9^(0) -9,(r) (g)9j(r)||oo < Cr. (31) 

It follows from (29), (30), and (31 )that € L^(i?p(xt); 

Now we compute 

^ 5]{c'(r)/'(i? + C(r)) : 9.(0) 0 9,(0) 

+ [/(^^ + C(?')) : 9'i(0) ® 9,(0) - f{^) ■■ gr{r)®gj{r)] 

• [g'i {r) ® 9, (r) + g^ (r) <S> g'j (r)] 

The analysis of the term l3g/r is similar to the one performed in (29), and thus 
d(3g/dr G L^(Bp(xt);M.^^^) follows from the boundedness of g[ and g '2 and (28a). □ 


3 . 4 . Simplified Expressions for the Interaction Integral Functionals 

We describe three pairs of material variation fields Sj and auxiliary fields /3^ = and 

for each pair we provide the simplified expressions of the interaction integral functional 
X[f3, (3‘^'^^,6'y] that results from substituting the two fields. In this section we have removed 
subscripts I, II from the auxiliary fields, as the following results are independent of the choice 
of the mode of interest, and doing so clarifies the presentation. 

We begin by stating two results used in obtaining the simplified expressions: (1) for traction- 
free auxiliary stress fields cr(/3®'™), such as cr{f3^^), and tangential material variation fields, 
such as we have 

J^TAN . - ^ -d^TAN . ^TF T (33) 

and (2) for compatible and divergence-free auxiliary fields, such as we have 

A(/3,/3°^°) (33) 


Justification (Equations (32) and (33)). We begin with (32). Recalling the expression (15) for 
r {f3, (3^'^^) we have, over 

Sj ■ T {(3, f3^'^^) = w i/3, (3^'^^) d 7 • n - d 7 • /3 ^o- (/3^"’^) n-S-f ^t. (34) 

Since we assumed that <57 is a tangential material variation field (i57 • n = 0) and because 
^(^aux) jg traction free (cr(/3®'™)n = 0 on '^^) then (32) holds. 

Next, we look at (33). Recall that, from (16), 

A(/3, 13^'^^) = (3 : V(Til3^'^^) - cr (/3) : (V/3“’^ )^ - /3 ^V • (Til3^'^^) +13^'^^ ^b. 
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Since we assumed that is divergence-free, the third term above vanishes. Furthermore, 

since we assumed that is compatible, there exists Bp(xt) \'rf—> such that 

^aux _ Exploiting the major and minor symmetries of the constitutive tensor C, we 
have 

/3 : Vcr (/3“’') - cr(/3) : ^ = /3 : C : VV# - (3 : C : VV$ = 0. 

Thus, (33) holds. □ 

We now present the simplified expressions for the functional obtained for each pair, as well 
as for the particular case of rectilinear cracks, in order to re-connect these results with what 
is commonly found in the literature. 


(1) Unidirectional material variation with divergence-free auxiliary fields. We set 6'y = 

and /3®'™ = Then (20) implies that the domain of integration of X 2 reduces to 

[Bp{xt) \ Bpj{xt)] \ ^. Substituting (33) in (17) then simplifies to 


X [/3, , 57 Uni] = f j^uNi . ^ ^DFC^ dS- [ (3°^^ ■ ^ 7 ^^! dV. 

-L 




o(xt)\'^ 

S(/3,/3°^®) : V57™i dV. 


(35) 


(2) Tangential material variation with divergence-free auxiliary fields. A slight variation of the 
previous pairing is the combination ^7 = ( 57 ^'’^'^ and ^9®'™ = Applying (34), (22), 

and (33) yields 


I [j3, ( 57 "^^^] = - [ [( 57 "^^^ • (3 ^cr n + < 57 ^^^ • ^t] dS 

[S {(3, ■ <57TAN] dV. 




(36) 


(3) 


Tangential material variation with traction-free auxiliary fields. Here we employ 6'y = 
( 57 TAN Applying (34), (22), (26) leads to 


I [f3, /3^^, ^7^^'"] = - / <57^^^" • f3'^^ ^ 


tdS 


f [S (/3, /3^^) : V,57T^n ^ ^tf^ . ^^tan j 

J Bp{xt)\V 


(37) 


(4) Locally rectilinear cracks. Finally it is worth noting that in the particular case of a locally 
linear crack geometry, i.e. r"(r) = 0, Vr € [0, p], 6'y = ( 57 '^an _ ^^uni _ ^dfc _ 

(3"^^, the interaction integrals of (35), (36) and (37) all simplify to 

I[f3, f3^'^^, Sy]= - [ Sy (3^'^^ ^tdS- [ "E{f3, f3^'^^) : VSy dV 

- [ (3^'^^^b ■ 5y dV 

J Bp{xt)\'^ 

which is the traditional expression of the interaction integral for a straight crack first 
introduced in [15] and commonly found in the literature. 

The presence of singularities in some of the factors in each one of the terms in (17) raises the 
question of whether the integrals therein are well-defined. For the choices of material variation 
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and auxiliary fields above, the three terms are in fact integrable. It is straightforward to see 
that I 2 is integrable, and the integrability of X\ and Xa is discussed below. 

Justification (Integrability of Xi ). We show the integrability of Ii by considering each term 
of T (cf. 15). For the first term of r, notice that w{(3, 0{r~^) and • n= 0{r) 

as r —)■ 0, and hence the first term in S'f -t asymptotically behaves as a constant near the 
crack tip. For the second term of r, we only need to consider the case Notice 

that on the crack faces, i? = ±7r — C{r), and therefore cos(i?/2), cos(3i?/2) = 0(r). Using the 
expressions for implies that 

: gfiO)(g)gfiO),(T{(3^^^): SaW ® 02 ( 0 ) = 0{r) 

on the crack faces near the crack tip. Moreover, since n = g2{r) — 92(0) + 0 (’") gi( 0 ) • 

92(0) = Oj have = 0(r). Finally, since /3 = close to the crack tip, we 

can conclude that /3^cr(/3'^^®)n ~ as r —>■ 0, and hence it is integrable. For the third term 

in T, if t e then Sj ■ /3®'™ as r —> 0, which is integrable as well. □ 


Justification (Integrability of I 3 ). As discussed in §3.4, we know A(/3,/3^^^) = ^b. If 

b G then (3^^^ as r —> 0. Therefore, for (3’^'^^ = X 3 is integrable. 

For = (3^^, we begin by taking advantage of (27) and the linearity of A in the second 
argument to write 

A {(3, (3^^) = A (/3, (3^^^) + A {(3, (3s) ■ 


But from earlier discussion about the regularity of (3^^, '^(^s = 0{r ^/^) as r —)• 0. Thus, it 
is straightforward to show that A {(3, (3s) = 0{r~^), and hence X 3 is integrable. □ 


3.5. Choosing the Interaction Integral Functional to Use 

Before introducing the numerical approximation of the above integrals it is worth making some 
remarks on which functional is best suited for a specific application. 

When the crack faces are loaded, a boundary integral over the faces has to be carried out 
irrespective of the auxiliary fields. For this particular problem it may be appealing to choose 
a pairing with such as (35) and (36). Doing so reduces the numerical complexity of the 

interaction integral as A greatly simplifies (and vanishes identically in the absence of body 
forces). 

If the crack faces are traction free, it can be appealing to compute the value of the interaction 
integral merely as a domain integral, as in (37). This eliminates the need to construct 
quadrature rules over the crack faces. Furthermore in the presence of body forces, the integrand 
A will be non-zero even with (3^^^, thus requiring the computation of the domain integral. 
For this particular case using (3"^^ will result in a computationally more efficient technique. 

Remark (Omission of unidirectional material variation with traction-free auxiliary fields). The 
pairing with (3'^^ is omitted because it provides no advantage over other pairings. In 

fact, because of 5^^^^, we have to compute the boundary integral Xi regardless of the loads on 
the body. Similarly, because of (3'^^, we have to perform the domain integral associated with 
the divergence of the reciprocal energy momentum tensor regardless of the loads on the body. It 
is thus apparent that for this particular pairing we do not eliminate neither Ii norl^, unlike 
for other pairings (when traction and body forces are zero). Therefore this pairing would result 
in computationally inefficient formulation, with no apparent advantage over other pairings. 


4. NUMERICAL COMPUTATION OF THE INTERACTION INTEGRAL 

In this section we are concerned with the computation of the interaction integral between 
any of the auxiliary fields and the solution u of Problem (I). The solution u and its gradient 
jd = Vu are going to be approximated by a convergent sequence of displacement fields {u^jh 
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and strain fields {j3^ = Vu^}hj respectively, or discrete solutions. We first give some general 
considerations on the expected conditions for convergence of the computed stress intensity 
factors that are independent of the method adopted to compute (3^. Then we particularize 
some of these results to a {/3^} that stems from a sequence of finite element approximations. 
Additionally, we discuss some minor changes needed when the exact domain needs to be 
approximated as well because of the presence of curvilinear cracks. The result of this section 
is an algorithm to compute T, summarized at the end of this section for readers interested in 
its implementation. 


4 . 1 . Approximation of the Interaction Integral 

Given a sequence of discrete solutions /3^ —)■ /3 in a sense to be specified later, it defines 
a sequence of values for the interaction integral ^ 7 ], and hence a sequence of 

approximate stress intensity factors 


^iji 






(38) 


For the approximate stress intensity factors to converge to the exact ones as h \ 0, 

it is enough for I to be continuous with respect to its first argument in the topology in which 
(3^ converges to /3. It is simple to see then that for the stress intensity factors computed with 
(37), it is enough to have —)■ /3 in because these functionals do not involve 

integration of /3 over In contrast, for the stress intensity factors computed with (35) and 
(36), we additionally need to request /3^ —> /3 in L^(^±; 


4 . 1 . 1 . Finite-Element-Based Approximations For sequences {u^}h constructed with some 
finite element spaces there is an important advantage of having a functional I continuous in 
its first argument. That is, the order of convergence of the stress intensity factors doubles the 
order of convergence of /3^ to f3 [32, 33, contain related results, and see §B], so the values of 
the stress intensity factors are a lot more accurate than the discrete solution itself. It is not 
difficult to check that I in (37) is continuous in its first argument in Therefore, 

we can conclude that if |j/3^ —)■ / 3 ||i 2 (e.^.R 2 x 2 ) < then j.j(/3^) — Kiji{l3)\ < for 

some C > 0, fc € N independent of h. 

The functional I given by (35) or (36) is not continuous in its first argument in 
because of the boundary integrals. As described in §B, the result that states that the order 
of convergence of jj should double that of (3^ does not apply in this case. Nevertheless, 
as shown later in the numerical examples, the rates of convergence seem to double as well for 
these two functionals. 

The values of k of the numerical methods used for the numerical examples in §5 are 0.5 and 
1, and thus these methods converge at the rates of 0{h) and 0{h?), respectively. In order to 
achieve higher order of accuracy within the context of finite element methods, it is necessary 
to make use of alternative methods that can accurately resolve the stress singularity, such as 
[I, 2, 3]. Furthermore, for curvilinear cracks, high-order approximations of the crack faces are 
needed to attain a corresponding order of accuracy of the method. 


4 . 2 . Discrete Interaction Integral Functional 

One of the delicate issues to be addressed in this section is the fact that for curvilinear cracks 
each discrete solution is computed on an approximation of the exact domain. The precise steps 
to handle the difference between exact and approximate domains in finite element methods are 
fairly standard, and hence are often skipped in the description of new methods. We decided 
to discuss this part with some additional detail here because of the presence of the boundary 
integrals. The uninterested reader could simply skip to the next section. 

For each h, the discrete solution is computed on a domain with crack faces We 
assume that as h \ 0 the approximate domain and the approximate crack faces and their 
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normal vectors converge to the exact ones^. For example, a standard isoparametric mapping 
will suffice. 

Because the discrete solutions are defined over different domains, the interaction integral 
functional T needs to be approximated as well (the integrals over and could not 
be computed for the discrete solutions otherwise). Thus, for each h we construct a discrete 
interaction integral functional At —>■ M, where is defined analogously to 

but considering as the domain of the problem. Then, given a sequence of solutions 
converging to the exact solution u in we expect , (^ 7 ] = X[/3'^, ( 57 ], 

for any of the and any ^7 C At, where f3^ = and /3® = Vm. Equivalently, 

letting the approximate stress intensity factors iCj jj : —> M be 




[/3\/3?"A,<57] 


(39) 


we expect lim/jx^oand limhx^o Kj = Kiji{(3‘^). These ideas are compactly 

shown in the following commutative diagram: 


, -rh 

/ 3 "-^ K^.ii 


\ 0 


h\0 


(3^ 


~^Ki, 


II 


The functional is defined as 

X^ [/3^/3-^^7] = E ^7•r^(/3^/3“")L^^ag 

- ^ [S (/S'*, : V57 + A (/j'*, 

g&G 

where Q and Q<^ denote the set of quadrature points over Bl^ and n Bp{xt), respectively. 
Each integration point 5 in ^ or Q<^ has position Xg and integration weight Wg. We assumed 
that all quadrature points over B^ belong to Bl^ n B<^ which is true for a small enough mesh 
size, to be able to evaluate (3^'^^, which is defined over B‘^. Additionally, we defined as an 
approximation to r given by 

-h ^aux^ ^ . ^(^aux o p)n O p - /3^ ^Cr(/3“" O p)n O p - ^t) O p, 

where p : 'rf^ 1 —>■ is the constant-radius projection of a point onto the crack: 

p(x) :=r(|x-a7t|). (41) 

This projection is well defined when and are close enough. Other projections are possible 
as well. This one is convenient, since it is also involved in the definition of 0^^. 

Remark (Appearance of p in the boundary integral). The functional ( 4 O) is an approximation 
to integrals over Bl^ and and the quadrature points of Q and Q<g belong to them. 


possible condition is that for each h there exists a one-to-one map G that 

converges to the identity in at a suitable rate, with det > e for some e > 0 uniformly in h, and 

such that B^ = and = This is a type of condition for finite element approximations, 

and it is simply a condition on the way the approximate domains are to be constructed; we will not need to 
explicitly construct to compute the interaction integrals. 

§For example, u — o 0 in 


<57] I 


(40) 
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The material variation and auxiliary fields are constructed over the exact domain and 
could be tangent or traction free to its boundary ^±, respectively, but not necessarily to its 
approximation Furthermore, the traction i is prescribed and only known over the exact 
crack '^±. To address this difficulty the auxiliary traction cr{l3°''^^)n, the auxiliary fields 
themselves, and the applied traction t are evaluated at their constant radius projection onto 
the curved crack p(a:) for x G H Bp{xf). Note as well that (^7 o p = ^7, since S'y depends 
only on r. 

Figure 6 shows an example for piecewise linear interpolations of the exact geometry along 
with its discrete approximation and the mapping p in (41). As the mesh is refined, the difference 
between p(a;) and x should go to zero, and the Jacobian of the mapping p : —>■ 'i^± should be 

very close to unity, thus permitting the composition in the boundary integral without introducing 
significant errors. 



Figure 6. Example of constant radius projection for piecewise linear interpolations of '^±. 


Remark (Convergence of the singular boundary integral). Recall that if the applied crack-face 
traction is bounded at the origin, we expect the boundary integral of (8) to possess a singularity 
at the crack tip as [3°''^^ oc for r —>■ 0. Integrating a singular function using standard 

Gaussian quadrature over a successively refined discretization was observed experimentally to 
lead to errors of the order ©(/i^/^) (see [2, Appendix B]). Therefore, in the particular case in 
which t is bounded and non-zero at Xt, it is necessary to address the numerical integration of 
the singular function in order to preserve the expected convergence rate. Here we computed the 
■singular integral 


[ 5-f ■ dS = r (,57 • (r)|r'(r)| dr 

Jo 

by pulling back the integrand (^7 •(r)|r'(r)| from [ 0 ,/?] to s“^([ 0 ,p]) through the 
map s(f) = (r^/p — f) q(f) + f with q of (18). Then we simply use the quadrature rule Q<^ over 
s“^([0,p]). Namely, if we let rg := |a;g — xf for all g G Gv, we compute the above integral as 





,57 

g&Q’g 


|s'(rg)|wg. 

sCs) 


The mapping s effectively performs a local change of variable r ^ r"^ which removes the \Jr 
singularity of the integrand thus allowing to recover optimal rates of convergence. The scaling 
of the 1/p ins serves to ensure that the mapping s is injective over [0,/)]. 
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4-3. Summary of the method 

We provide here a very concise summary of the method for the reader seeking a guideline for 
a rapid implementation. 

The calculation of the stress intensity factors can be summarized in the following steps: 
Step 1: Compute the approximation (3^ to the gradient of the solution of (1). 

Step 2: Construct 5^ to be either or < 57 '^^^ from (19) or (21) respectively. 

Step 3: Construct to be either /Sjj/ or /3/from (25) or (24) respectively. 

Step 4- With the pair ^ 7 , as well as b, t and (3^ use (40) to compute the value of . 

Step 5: Compute the value of jj with the above following (38) (or (39) if appropriate). 

We recapitulate int Table I the simplifications of each integrand associated with each choice 
of pairing of and ^ 7 . 

Table I. Recapitulation of simplifications associated with the choice of fields. Omission of terms (-) 

stands for no simplification. 


Fields 


I 2 {B^) 

^3 {B^) 

^7Uni^^dfc 

- 

= 0, Vr < Pi 

(57UNI ./3DFC 

^7TAn^^dfc 

S-y ■ [f3'^^ - (3'^ ^(T(/ 3 °^®)n] 

- 

^7TAN . ^DFC T5 

^7TAn^^tf 

5-f ■ 

- 

- 


5. NUMERICAL EXAMPLES 

We next verify the proposed method through two examples. For each we provide comparisons 
with analytical solutions. The first problem is concerned with a circular arc crack in an infinite 
medium subjected to far-held stresses. The second problem involves a power function crack 
subjected to crack face tractions and body forces. 

For each example we compare the convergence of the stress intensity factors for lower 
order methods, namely traditional continuous Galerkin hnite element methods for piecewise 
polynomial shape functions fc = l,2, and for the higher order discontinuous Galerkin 
extended hnite element method (DG-XFEM) [2]. Both methods are recapitulated in §5.1. 

As discussed in §4, the interaction integral, and hence the stress intensity factors, are 
expected to converge at twice the rate of the derivatives of the solution. Thus we are expecting 
to observe convergence of the order 0 {h) for lower order methods (whose derivatives converge 
as 0{h^'^)) and 0{h3) for the higher-order DG-XFEM method (whose derivatives converge as 
0 (h^)), where h is the maximum diameter of a triangle in each mesh in the family of meshes 
under consideration. 

In the following examples we will provide systematic convergence curves of the error in the 
solution and in the computation of the stress intensity factors. Tabulated errors and computed 
convergence rates will accompany the above. 

We will present two error measures of the solution, one over the interior of the domain and 
the other over the crack faces. The error in the solutions over the interior of the domain will be 
measured as the L^-norm of the error in the gradient of displacements over and that over 
the crack faces will be measured as the L^-norm of the error in the gradient of displacement 
weighted by r over Namely, with /3'^ denoting the analytical solution of the gradient of 
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displacement fields, we will consider as error measures 


W - 




J 


1/2 


(42) 


and 


/3-rop 


L^{V^,r) 


■ n 1/2 

E,, - Pf, op)\ds =11 (/3'* - r o p) . 


(43) 

In (43) the analytical gradient of displacements are evaluated at their constant radius 
projection onto the exact geometry as discussed earlier in §4.2. 


Remark (Appearance of ^/r in the L^('^jl) norm). The appearance of the ^/r factor is related 
to the scaling of the factors that multiply (3^ in the integrand of Ti. Namely (3^ appears as 
/3^ ^(T(/3“™)n and (3^ : (T{f3°‘'^^)6'y ■ n. In the former we have cr{f3‘^^^)n ~ r as r —)■ 0, as 
previously discussed in section %3.4- In the latter we need to consider the scaling of Sj ■ n, 
which is either (57™^ • n = 0 or ■ n ^ r, as well as the scaling cr(/3““'^) ~ Ij\/r, as 

r —> 0. Hence f3^ in the latter case is multiplied by a factor that scales as \/r as r —>■ 0. Thus, 
only the rate of convergence of ^/rf3^ is needed to evaluate the rate of convergence of Ti . 


The error in the stress intensity factors will be measured by the normalized absolute value of 
the error in the computed stress intensity factors. Namely, let Kjjj := Kjjj[(3^] be computed 
with (38) (or (39)) and Kf be the exact (analytical) stress intensity factors. We will be 
concerned with the behavior of 

I ffh _ T^e I 
I /d® I 

We will also present for each example the value of the computed stress intensity factors for 
various values of p, that is, for different supports for 5^. As the interaction integral in (8) 
is independent of p, we would like to test the independence of the computed stress intensity 
factors on the support of 5'^. 

Lastly we remark that for each example we set the material constants to A = 277.77, /r = 2500 
{E = 1000, n = 0.2) and we assumed a plane strain state. 

5.1. Numerical Solution of the Elasticity Problem 

We consider two types of finite element methods over a family of meshes of triangles. In the 
following, the superscript (•)^ will denote quantities associated with the discrete approximation 
of the problem. For each mesh in the family, the domain B is approximated by 
the collection of open, straight triangles T®. Let V denote the set of all vertices in the mesh. 
Each mesh in the family conforms to the crack, namely, a node sits at the crack tip, and there 
is no edge with its two vertices on different sides of the crack. To handle the displacement 
discontinuity across the crack, vertices that lie on are duplicated, and so are edges whose 
two vertices lie on . The union of these edges on either side of the crack forms the piecewise 
linear approximation to to±, and we set Bl^ = For convenience, we define 

= {a e V|a;, € daB}, 

where Xa represents the position vector of vertex a. In the following examples, we let 
Na&H^{Bl^) be the shape function associated with node a G V, A: = I,2, such that 
Na{xb) = dab for all a,b gY, where 6 ab is the Kronecker delta. Of course, for the piecewise 
quadratic case k = 2, mid-edge nodes are added . 

The two methods adopted here are: 
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(1) Standard finite element method. We seek an approximate solution G S^, with 
= {EaeV^aWa = u’^ G {B^;R^)\u\xa) = u{xa),ya G V'"} . 
We further let 

v'* = {Eaev^a<5Ma = Su’^ G {B!^-,R^)\ 6 u\xa) = 0,Va e V'^} . 
The numerical approximation of u is obtained by hnding G such that 



: V5u^ dV 


b ■ Su^ dV + 




Su^ dS, 


where 

=S/U^ = '^Ua®VNa. 
oev 


(2) Discontinuous-Galerkin extended finite element method. Here we recapitulate the method 
proposed in [2] with slight improvements. Let h < (l/2)[dist(a7t, i9H) — p] and be such 
that 

p h < rc < dist(£C(, 9H) — h, 

and _ 

B^ = [j{T-\aye4T-nBrAxt)]>0}, B'i =W\^, 

be the enriched and unenriched regions, respectively. Then we set 

= {aGN\xaGBl), = {aGY\xaGBj^}. 

Hence, there are nodes that belong to both Y^ and Y^. In fact, let Tf = dBj^, then 


V®nv'^ = {aGV|£Caerf}. 


The discontinuous Galerkin extended finite element method (DG-XFEM) is built on the 
following set: 

S’^ = {u’^ G {B%- m2) \u^ = kiu^ + kuu^^ + Eaev- in , ku ku & K; 

vfi = EaeVV in U^{Xa) = u{Xa),ya G . 

The corresponding test space is given by: 

= {Su^ G {B^-, m2) \6 u'^ = Skiu^ + 6kiiu^^ + ^o.5u^ in , 6 ki, 5kn G M; 

5vfi = J2aeVo in Bj^,6u^{xa) = 0,Va G V'^} . 

Therefore, the kinematics of a typical function G is independent in Bj^ and ; a 
discontinuity across Vf arises which is defined as 

|m ] = M IgE ~ Igv = [klU + kljU + EaeV^nVV {u^ — 'It^)]pE ■ 

LhJI^LhJT^ h 

This discontinuous is handled through a DG-derivative Ddg ■ —> YhS^ + 

DoG'-u^^Ynu'^ + Rilu^l), 
where YhU^ = Yu^ in each T®, i?(|M^]) is such that 



: dS, Vic'* e 
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where 




and on rf 


{-n = d 


w 


'’/i J r 


w 


The solution to the problem stated in §2.1 is approximated by: Find such that 


/ cr (/3'') : DdgSu^ dV + 2na [ R\ 

JBh JBh 


M dV 


b ■ (5m'* dV + 


Id-rBU^^ 


t ■ 6u^ dS, ySu'^ e V'*, 


where a can be any positive real number, and 

/s'* = Ddgu^ 


We conclude the section by remarking that the approximate domain of integration of the 
interaction integral for the particular choice of the method is given by the subset of elements 
with at least one vertex that lies within Bp{xt) which we denote by K. Refer to Fig. 7 for an 
illustration of the above. 



Figure 7. Example of a subset of elements IK in a finite element mesh. The elements in the shaded 
region are the elements over which the interaction integral is computed. 


Furthermore, we exploit the quadrature rule constructed over each element Qt<^ and its 
boundary Qqt’^ to form Q and respectively. Namely we let the numerical interaction 
integral, in this specific setting of finite element methods, become 

T=GK,|0T=n<^^|>0 g^GdT^ 

-EE [S(/3\/3“") : V57 + A(/3\/3"'*") -57] \^Wg. 

g^Q'j'e 


Prepared using nmeauth.cis 






COMPUTING STRESS INTENSITY EACTORS EOR CURVILINEAR CRACKS 


25 


5.2. Circular Arc Crack 

We consider an infinite plate with a circular arc shaped crack subjected to uniform tension 
from infinity. The analytical solution was derived in [34], and a recapitulation of the solution 
can be found in [35]. The resulting stress intensity factors for uniform far field tension loading 
are (see, e.g., [36]) 



where R is the radius of the circular arc crack, a is half the angle subdued by the crack, and 
a is the far field tension as shown in Fig. 8. 

Only a finite subdomain was considered and exact tractions were specified on the boundaries. 
Given the symmetry of the problem, only half of the subdomain was modeled and appropriate 
symmetry boundary conditions on the axis of symmetry were specified. Figure 9 shows a 
representation of the modeled subdomain and boundary conditions. For the simulations we 
took a = 7r/2, i? = 1, the modeled domain was given hy B = [0,2] x [—1.5,0.5], and the crack 
centered at the origin. 

To establish the accuracy of the methods, the solution was computed for different levels 
of refinement of the discretized domain. The meshes were generated by conforming recursive 
subdivisions of the coarsest mesh to the exact geometry. 

The error measures (42) and (43) were observed to decrease as for the lower-order 

method, and as 0{h) for the second-order method. Figure 10 shows the convergence plot of 
the solution, and Table 11 summarizes the error as well as the computed rates of convergence. 

As expected, the error in the stress intensity factors are observed to converge with order 
0{h^) and 0{h?) for the lower- and higher-order methods, respectively. Figure 11 provides the 
convergence curves for the stress intensity factors using the three pairings of material variation 
and auxiliary fields. Errors and computed rates are reported in Table Ill. 

Lastly we show that the evaluation of the interaction integral is independent of the support of 
5^. To this end, Fig. 12 shows the error in computed stress intensity factors of the most refined 
mesh for five values of p/pmaxj ranging from ~ 0.7 to 1 with pmax = 0.5. The independence of 
the interaction integral on the choice of the support of the material variation field is apparent 
from these results. 





>-oo 


= crl 


q 

CtJ 

II 



Figure 8. The circular arc crack problem. 


Figure 9. Modeled subdomain. Here is the exact 
traction on a face. 
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Table II. Convergence rates of the derivatives of the solution for the circular arc crack problem 


(a) Domain convergence 




pi 


p2 


DG - XFEM 


hmax^/ h 

Err. 

0 

Err. 

0 

Err. 

0 

1 

0.00055 

- 

0.00025 

- 

0.00028 

- 

2 

0.00037 

0.57 

0.00019 

0.43 

0.00015 

0.92 

4 

0.00026 

0.52 

0.00013 

0.51 

0.00007 

0.97 

8 

0.00018 

0.51 

0.00009 

0.51 

0.00004 

0.96 

16 

0.00013 

0.51 

0.00006 

0.51 

0.00002 

0.98 


(b) Trace convergence 


11/3^ — /3® o 



pi 


p2 


DG - XFEM 

^max/ h 

Err. 

0 

Err. 

0 

Err. 

0 

1 

0.00055 

- 

0.00025 

- 

0.00028 

- 

2 

0.00037 

0.57 

0.00019 

0.43 

0.00015 

0.92 

4 

0.00026 

0.52 

0.00013 

0.51 

0.00007 

0.97 

8 

0.00018 

0.51 

0.00009 

0.51 

0.00004 

0.96 

16 

0.00013 

0.51 

0.00006 

0.51 

0.00002 

0.98 



Figure 10. Convergence of the solntion. 


5.3. Power Function Crack 

The second example we consider is the one of the power function crack '^ = {{x, a:^)| x € [0,1]} 
loaded by a force field b and crack face traction t, see Fig. 13. 

The exact stress field is constructed by a superposition of a singular stress field it® with a 
bounded field 0-5 as 


cr 


e 


= d-® + 0-5. 
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h/hfjiax 

(c) (^"■“’',57) = (/3pFC,5-),UNi)^ 




h'l hjjifix. 






(e) (^=■“’',57) = (/3pFC,5-),TAN)^ (f) ^ (^DFC^^^TAN)^ 


Figure 11. Convergence of the stress intensity factors for the circular arc crack 
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Table III. Convergence rates for stress intensity factors 


(a) 


Traction free auxiliary fields = 


/3^^) and tangential material variation {S~f = 


pi 

p2 

DG- 

XFEM 

Ki 

1 

Kii 

Ki 

1 Kii 

Ki 

Ku 


^max./ h 

Err. 

O 

Err. 

O 

Err. 

O 

Err. 

O 

Err. 

O 

Err. 

O 

1/1 

5e-02 

- 

5e-02 

- 

4e-03 

- 

2 e-02 

- 

5e-03 

- 

2 e-02 

- 

1/2 

2 e-02 

1.75 

3e-02 

0.82 

5e-03 

0.11 

9e-03 

1.24 

9e-04 

2.45 

5e-03 

1.99 

1/4 

8e-03 

1.00 

le-02 

1.03 

2e-03 

1.08 

4e-03 

1.06 

3e-04 

1.69 

le-03 

1.91 

1/8 

4e-03 

1.08 

7e-03 

1.02 

9e-04 

1.19 

2e-03 

0.94 

le-04 

1.48 

3e-04 

1.91 

1/16 

2e-03 

0.99 

4e-03 

1.03 

5e-04 

0.92 

le-03 

0.92 

2e-05 

2.20 

7e-05 

2.23 


(b) Divergence-free and unidirectional material variation (S'y = ). 


pi 

p2 

DG- 

XFEM 

Ki 

1 

Ku 

Ki 

1 

Ku 

Ki 

Ku 


^max/ ^ 

Err. 

O 

Err. 

O 

Err. 

O 

Err. 

O 

Err. 

O 

Err. 

O 

1/1 

le-02 

- 

3e-02 

- 

7e-03 

- 

le-02 

- 

7e-03 

- 

le-02 

- 

1/2 

2e-03 

2.54 

2 e-02 

0.44 

9e-04 

2.93 

5e-03 

1.23 

le-03 

2.59 

le-03 

3.04 

1/4 

le-03 

0.91 

9e-03 

1.05 

5e-04 

0.92 

2e-03 

1.33 

le-04 

3.37 

2e-04 

2.42 

1/8 

7e-04 

0.76 

5e-03 

0.99 

3e-04 

0.78 

le-03 

0.81 

6e-05 

0.92 

8e-05 

1.52 

1/16 

3e-04 

0.99 

2e-03 

1.00 

9e-05 

1.57 

5e-04 

1.05 

8e-06 

2.77 

3e-05 

1.76 


(c) Divergence-free and tangential material variation {5~f = ). 




p2 


DG - XFEM 


Ki 


Kii 


Ki 


Ku 


Ki 


Ku 


^max/ h 

Err. 

O 

Err. 

O 

Err. 

O 

Err. 

O 

Err. 

O 

Err. 

O 

1/1 

4e-02 

- 

le-02 

- 

6e-03 

- 

7e-03 

- 

6e-03 

- 

8e-03 

- 

1/2 

le-02 

1.74 

le-02 

0.20 

4e-03 

0.56 

3e-03 

1.41 

8e-04 

2.85 

2e-03 

2.43 

1/4 

7e-03 

1.00 

5e-03 

1.04 

2e-03 

1.06 

le-03 

0.94 

2e-04 

2.09 

4e-04 

1.88 

1/8 

3e-03 

1.08 

3e-03 

0.90 

le-03 

0.82 

6e-04 

1.31 

5e-05 

2.11 

9e-05 

2.23 

1/16 

2e-03 

0.97 

le-03 

1.03 

5e-04 

1.13 

2e-04 

1.23 

le-05 

2.07 

2e-05 

2.14 


(M 



O' 

W 


10- 


10 - 


10“ 


0.6 


< 1_<1 (57 = g{r)gi(0) & divCT(/3S) = 0 
57 = q{r)gi{r) k divCT(/3f/j) = 0 
0-0 ^7 = q{r)gi{r) k (T{/3f‘fj)n = 0 


< 3 - 

[Ib 


0.7 


-O 


-<b 

-o- 


-o- 


0.8 


0.9 


PjPn 


-□ 


1.0 


1.1 


Figure 12. Contour independence of the interaction integral 
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The field is constructed as 


= Ki + O'// )9*(0) ® 9i(0), 

where are given in §A and are evaluated for values of -d C [— tt — C(r), tt — (^(r)], as 

discussed in §2.2. The bounded stress field 0-5 is constructed as 

cTf, = xe^ ®ex+ ysy 0 By. 


Note that 

V • = V • d-® + V • (Tb = V • (Tb = -b, 

where 

b = —{e^ + By). 

It is worth remarking that the stress intensity factors of cr® will correspond to those of the 
singular field d"® without any perturbation from the bounded field. In fact, given that both 
limr_i.o 1 linir-s-o \pr<^b — 0 exist, we have 

lim y/rcr^^ = lim + CTb) = lim y/rcri // + lim y/rcri, = lim \/rcrj ij. 

r —>-0 1 —^0 r —>-0 ’ r —>-0 r —>-0 

For our example we take the stress intensity factors of d"® to be Kf = Kfj = I. 

Figures 13 and 14 show the schematic of the problem and the modeled domain with the 
applied boundary conditions, respectively. For the simulations the modeled domain was given 
by ,8= [0,2] X [-0.25,1.75]. 

Like for the previous example, we computed the solution for several levels of refinement and 
investigated the convergence of the computed stress intensity factors. 

The error measures (42) and (43) were observed to decrease as and 0{h) for the 

first- and second-order methods, respectively. The values are plotted in Fig. 15 and the errors 
and the computed rates of convergence are tabulated in Table IV. 

The stress intensity factors were observed to converge to the analytical value as 0{h) and 
0{h^) when using the solution of the first- and second-order methods, respectively. The values 
of the error in the stress intensity factors are plotted in Fig. 16, and the errors, as well as the 
computed convergence rates, are provided in Table V. 

Lastly, in Fig. 17 we illustrate the independence of the interaction integral on the size of the 
support of Sj by plotting the stress intensity factors for five values of pjpmax ranging from 
~ 0.7 to 1, for Pmax = 0.5. 
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t = (T^n 

Figure 13. The power function crack problem 
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Figure 14. Modeled subdomain 


Table IV. Convergence rates of the derivatives of the solution for the circular arc crack problem 


(a) Domain convergence 




pi 


p2 


DG - XFEM 


^max/ h 

Err. 

0 

Err. 

0 

Err. 

0 

1 

0.00080 

- 

0.00040 

- 

0.00026 

- 

2 

0.00061 

0.41 

0.00032 

0.36 

0.00012 

1.09 

4 

0.00043 

0.50 

0.00022 

0.50 

0.00006 

1.03 

8 

0.00031 

0.47 

0.00016 

0.49 

0.00003 

1.00 

16 

0.00022 

0.50 

0.00011 

0.50 

0.00002 

1.00 


(b) Trace convergence 


11/3^ — /3® O p\\L^('^h^r) 



pi 


p2 


DG - XFEM 

^max/ h 

Err. 

0 

Err. 

0 

Err. 

0 

1 

0.00066 

- 

0.00037 

- 

0.00035 

- 

2 

0.00050 

0.40 

0.00029 

0.34 

0.00013 

1.39 

4 

0.00036 

0.48 

0.00021 

0.48 

0.00006 

1.06 

8 

0.00026 

0.43 

0.00015 

0.48 

0.00003 

0.97 

16 

0.00019 

0.49 

0.00011 

0.50 

0.00002 

0.99 
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(a) Convergence in the norm 


(b) Convergence in the norm 


Figure 15. Convergence of the solution. 
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Table V. Convergence rates for stress intensity factors 


(a) Traction free auxiliary fields 


/3^^) and tangential material variation ((5'7 = 


pi 

p2 

DG - XFEM 

Ki 

Kjj 

Ki 

Kii 

Ki 

Kii 


^max/ h 

Err. 

O 

Err. 

O 

Err. 

O 

Err. 
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Figure 16. Convergence of the stress intensity factors for the power function crack problem 
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Figure 17. Contour independence of the interaction integral 
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A. MODES I AND II ASYMPTOTIC SOLUTIONS 


We recall below the components of displacement, gradient of displacements and the stress fields 
for a straight crack lying on the axis i? = ±7r,r € [0,oo) as derived in [14]. The components 
are given for a set of right-handed orthonormal basis with the 1 axis aligned with the crack. 

For K = 3 — 4^ for plane strain and k = {3 — 4 ^)/(\ -\- v) for plane stress, the displacements 
are given by 

The gradient of displacements are given by 
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Lastly the stress components are given by 
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Note that the above fields satisfy 


Ki[Vu^] = l, Kii[Vu^]=0, 

Ki[Vu^^] = 0 , Kii[Vu^^] = 1 . 


B. CONVERGENCE OF A CONTINUOUS AFFINE FUNCTIONAL 

In this appendix we first state and prove a proposition about the convergence of linear and 
continuous functionals of a convergent family of solutions. This proof is essentially adapted 
from similar results in [33, 32]. We next apply this result to the interaction integral functional 
(37). 

Proposition. Let V he a Hilbert space and C V be its finite dimensional approximation. 
Let a:UxU —be a bilinear, continuous and coercive form with a{u,v) < Ci||m||\/||u||v 
for all u,v G V. Let F : U —)■ M 6 e linear and continuous functional, and let G ■. V ^ W he an 
affine and continuous functional. Take u to he the solution to a{u,v) = F{v), Vv G V and 
the solution to a{u^,v^) = F{v^), . Let w G V be the unique member ofV such that 

a{v,w) = G(v) — G(0), Vu G V. We further assume that there exist positive real numbers G 2 
and G 3 independent of h such that jju — < ^ 2 / 1 ^ and inf^hgyh \\w — w^jjy < G^h^^. Then 

there exists G independent of h such that 

|G(u)-G(m'‘)| <G/r2fc. 

Proof 

From the definition of w, 

G{u) — G [u^) = a{u — u^, w) . 

Furthermore note that for any G V^, 

a{u — u^, w) = a (u — w — + a (u — u^, w^) = a{u — , w — , 

where we have taken advantage of Galerkin orthogonality, i.e., 

a{u- u^, w^) = a {u, w^) - a (w\ w’^) = F {w^) - F {w’^) = 0. 

Therefore we have 

|G('u) — G (u^) I = \a{u — ,w — | < Gi ||m — Hw — . 

Since is arbitrary, we have 

|G(u)-G(u'‘)| < Gi inf |U - ^'*[1 < GiGaGg/i^^ 

Taking G = G 1 G 2 G 3 yields the conclusion. □ 

The application of this proposition to the interaction integral functionals here requires some 
additional work to account for the difference between domains in curvilinear cracks, and the 
use of quadrature rules. However, disregarding these differences, and assuming that 
and that exact quadrature is adopted, we have = I, and for the standard finite element 

method in §5.1, C \/h. Now, in its first argument, T[I3, (3^^, is affine 

and continuous in so we set G{u) = Z[Vm, (c.f. (37)) and can use 

the above proposition. Since for the standard finite element method it is known that there 
exists G > 0 independent of h such that [jw — < G/i*, V/i, then 

|G(t 6 )-G(M'‘)| <Gh 2 '= 
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for some C € M+. 

This proposition is not directly applicable to the discontinuous Galerkin method in §5.1, 
because in this case it is also necessary to account for the use of an approximation space that 
does not conform to H^. Finally, the two functionals G{u) = ^ 7 ] in (35) and (36) 

are not continuous in because of the evaluations of Vtt on the crack faces, so we 

cannot directly apply the above result. 
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